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ABSTRACT 

We discuss a numerical model for black hole growth and its associated feedback pro- 
cesses that for the first time allows cosmological simulations of structure formation to 
self-consistently follow the build up of the cosmic population of galaxies and active 
galactic nuclei. Our model assumes that seed black holes are present at early cosmic 
epochs at the centres of forming halos. We then track their growth from gas accretion 
and mergers with other black holes in the course of cosmic time. For black holes that 
are active, we distinguish between two distinct modes of feedback, depending on the 
black hole accretion rate itself. Black holes that accrete at high rates are assumed to 
be in a 'quasar regime', where we model their feedback by thermally coupling a small 
fraction of their bolometric luminosity to the surrounding gas. The quasar activity 
requires high densities of relatively cold gas around the black hole, as it is achieved 
through large-scale inflows triggered by galaxy mergers. For black holes with low accre- 
tion rates, we conjecture that most of their feedback occurs in mechanical form, where 
AGN-driven bubbles are injected into a gaseous environment. This regime of activity, 
which is subdominant in terms of total black hole mass growth, can be identified with 
radio galaxies in clusters of galaxies, and can suppress cluster cooling flows without the 
requirement of a triggering by mergers. Using our new model, we carry out TreeSPH 
cosmological simulations on the scales of individual galaxies to those of massive galaxy 
clusters, both for isolated systems and for cosmological boxes. We demonstrate that 
our model produces results for the black hole and stellar mass densities in broad 
agreement with observational constraints. We find that the black holes significantly 
influence the evolution of their host galaxies, changing their star formation history, 
their amount of cold gas, and their colours. Also, the properties of intracluster gas are 
affected strongly by the presence of massive black holes in the cores of galaxy clusters, 
leading to shallower metallicity and entropy profiles, and to a suppression of strong 
cooling flows. Our results support the notion that active galactic nuclei are a key in- 
gredient in cosmological structure formation. They lead to a self-regulated growth of 
black holes and bring the simulated properties of their host galaxies into much better 
agreement with observations. 

Key vifords: methods: numerical - black hole physics - galaxies: formation - galaxies: 
clusters: general - cosmology: theory 



1 INTRODUCTION 

It is now widely believed that most if not all galaxies with 
a spheroidal component harbour a supermassive black hole 
(BH) in their centres. Interestingly, the masses of these cen- 
tral BHs are found to be tightly linked with the stellar prop- 
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erties of their host galaxies, as expressed, e.g., in the correla- 
tion of bulge velocity dispersion with BH mass (Ferrarese & 
Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002), 
or in the relation of bulge stellar mass (Kormendy & Rich- 
stone, 1995; Magorrian et al., 1998; Marconi & Hunt, 2003; 
Haring & Rix, 2004) with BH mass. The existence of these 
relationships indicates that the formation and evolution of 
galaxies is fundamentally influenced by the presence of BHs, 
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and vice versa. We thus also expect that the environment 
and the cosmological evolution of galaxies will affect the way 
BHs grow. 

In fact, there is a plethora of observational and theoreti- 
cal studies that suggest that several different channels for in- 
teraction of BHs with their surroundings exist. At high red- 
shift, mergers of gas-rich galaxies occur frequently and fun- 
nel copious amounts of cold gas towards the central regions 
of galaxies, such that the embedded BHs can reach high gas 
accretion rates. The radiation energy associated with the 
accretion can support the enormous luminosities of pow- 
erful quasars. Theoretically it has been hypothesized (Silk 
& Rees, 1998; Fabian & Iwasawa, 1999; King, 2003) that 
quasars produce high velocity winds, which are expected 
to affect the properties of the host galaxy. The presence of 
quasar induced outflows has been observationally confirmed 
in a number of cases (e.g. Chartas et al., 2003; Crenshaw 
et al., 2003; Pounds et al., 2003), and has first been demon- 
strated in simulations of merging galaxy pairs by Di Mat- 
teo et al. (2005). Several numerical studios dealing with BH 
microphysics (e.g. Proga, 2003; McKinney, 2006) also pre- 
dict existence of quasar outflows. Moreover, it appears that 
tidally disrupted galaxies arc preferentially associated with 
AGN activity (for a review see Barnes & Hernquist, 1992). 
Quasar activity hence appears to be directly linked to merg- 
ers of galaxies, and should represent the dominant mode of 
mass growth in the BH population. 

Indeed, using semi-analytic models of galaxy formation, 
Kauffmann & Haehnelt (2000) have demonstrated that BH 
growth associated with mergers in CDM models can repro- 
duce many properties of the observed quasar population as 
well as the inferred BH mass density today (see also Volon- 
tcri et al., 2003). Based on the detailed hydrodynamical sim- 
ulations of BH growth in galaxy mergers (Di Matteo et al., 
2005; Springel et al., 2005a; Robertson et al., 2006a,b; Cox 
et al., 2006a,b) that have recently become available, Hopkins 
et al. (2005, 2006a,c) have proposed a comprehensive picture 
of a unified, merger-driven model for the origin of quasars 
and their relation to spheroid formation, which also implies 
the existence of a 'black hole fundamental plane' (Hopkins 
et al., 2007a). Also, rapid merging of the gas-rich progenitor 
systems of rare, massive galaxy clusters has been shown (Li 
et al., 2006) to be a viable formation path for supermassive 
BHs that are as massive at « ~ 6 as those seen in luminous 
high-redshift SDSS quasars (Fan et al., 2001). 

However, there also appears to exist another chaimcl of 
BH interaction with host galaxies, which is neither related 
to powerful quasar activity nor associated with galaxy merg- 
ers. Evidence for this interaction can be seen in a number 
of local elliptical galaxies and central cluster galaxies, which 
contain X-ray cavities filled with relativistic plasma (Birzan 
et al., 2004; McNamara et al., 2005; Forman et al., 2006; 
Fabian et al., 2006) while harbouring seemingly "dormant" 
BHs. These X-ray depressions, often referred to as 'bubbles', 
are thought to be inflated by relativistic jets launched from 
the central BH. Even though the radiative output from the 
central BH is not significant, the associated mechanical lu- 
minosity can be very important in these systems. There has 
been considerable effort from the theoretical point of view 
(e.g. Binney & Tabor, 1995; Churazov et al., 2001, 2002; 
Quilis et al., 2001; Ruszkowski & Begehnan, 2002; Dalla Vec- 
chia et al., 2004; Kawata & Gibson, 2005; Sijacki & Springel, 



2006; Thacker et al., 2006; Okamoto et al., 2007) to under- 
stand the relevance of this mode of AGN feedback, which is 
widely considered to be leading candidate for resolving the 
cooling flow problem in groups and clusters of galaxies. 

Given the complex physics of AGN-galaxy interactions, 
is it possible to construct a simple unified model that ac- 
counts for the different modes of BH feedback in a cosmolog- 
ical framework? First attempts in this direction have already 
been made (Churazov et al., 2005; Croton et al., 2006; Mer- 
loni & Heinz, 2006), motivated by the observational findings 
of X-ray binaries (Fender et al., 1999; Gallo et al., 2003). In 
particular, it has been shown that X-ray binaries switch be- 
tween two states: in the so-called "low/hard" state, a steady 
radio jet is present and the hard X-ray spectrum is observed, 
while in the "high/soft" state, the jet vanishes and the X- 
ray spectrum shows a soft, thermal component. The transi- 
tion between these two states is regulated by the accretion 
rate onto the BH itself, where the threshold value is of the 
order of 10"^ - 10"^AfEdd. The "high/soft" state can be ex- 
plained by the standard, optically thick and geometrically 
thin accretion disc (Shakura & Sunyaev, 1973), with BH ac- 
cretion occurring at high rates and in a radiatively efficient 
mode. The "low/hard" state on the other hand corresponds 
to optically thin, geometrically thick, and radiatively ineffi- 
cient accretion, as described by the theoretical ADAF and 
ADIOS solutions (Narayan & Yi, 1994; Blandford & Begel- 
man, 1999). 

The above suggests a rather simple, yet attractive sce- 
nario for distinguishing between different modes of BH feed- 
back in models for the cosmological evolution of active galac- 
tic nuclei (AGN): We shall assume that for high accretion 
rates, a 'quasar-like' feedback occurs, while for states of low 
accretion, mechanical bubble feedback applies. It is clear 
that the simplicity of this model will not allow it to ex- 
plain all kinds of AGN feedback phenomena, e.g. powerful 
radio galaxies that accrete at very high rates, as found in 
some proto-cluster environments, are not well represented 
in this simple scheme. Moreover, even though the physics 
of X-ray binaries is expected to be quite similar to the one 
of AGN (Heinz et al., 2005), an analogous transition be- 
tween "high/soft" and "low/hard" states in AGN is obser- 
vationally still not well established, although there are some 
encouraging observational findings (Maccarone et al., 2003; 
Kording et al., 2006) in this direction. 

With these caveats in mind, we explore in this study 
a new 'two-mode' AGN feedback scenario in fully self- 
consistent cosmological sirrmlations of structure formation. 
This complements and extends our study of cosmological 
simulations with quasar feedback in Di Matteo et al. (2007), 
where we did not account for the "radio" mode. Within our 
numerical model we represent BHs with coUisionless "sink" 
particles and we adopt subresolution methods to compute 
their gas accretion rate, estimated from a Bondi prescrip- 
tion. Also, we allow BHs to grow via mergers with other 
BHs that happen to be in their vicinity and that have suf- 
ficiently low relative speeds. We assume that BH feedback 
is composed of two modes, £is motived above, and we track 
their growth and feedback with cosmic time. We consider 
a vast range of objects harbouring BHs, from galaxies with 
~ IO^Mq stellar mass to massive galaxy clusters with to- 
tal mass ^ 10^^ Mq. Even though our BH growth model 
is rather simple and crude due to the inevitable numerical 
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limitations, it represents the first attempt to include AGN 
feedback effects in cosmological simulations of structure for- 
mation. As we shall see, this approach produces significant 
improvements in the properties of simulated galaxies. 

The outline of this paper is as follows. In Section 2, 
we describe our numerical method for incorporating the BH 
growth and feedback in cosmological hydrodynamical simu- 
lations. In Section 3, we test our model in simulations of iso- 
lated galaxy clusters, and explore the numerical parameter 
space. Wc then discuss cosmological simulations of galaxy 
cluster formation and evolution, subject to AGN feedback, 
in Section 4, while in Section 5 we consider cosmological sim- 
ulations of homogeneously sampled periodic boxes. Finally, 
we summarize and discuss our results in Section 7. 



2 METHODOLOGY 

We use a novel version of the parallel TreeSPH-code 
GADGET-2 (Springcl, 2005; Springcl ct al, 2001b) in this 
study, which employs an entropy-conserving formulation of 
smoothed particles hydrodynamics (Springel & Hernquist, 
2002). Besides following gravitational and non-radiative hy- 
drodynamical processes, the code includes a treatment of 
radiative cooling for a primordial mixture of hydrogen and 
helium, and heating by a spatially uniform, time-dependent 
UV background (as in Katz et al., 1996). Star formation 
and associated supernovae feedback processes axe calculated 
in terms of a subresolution multiphase model for the ISM 
(Springcl & Hernquist, 2003a). Wc use a simple prescription 
for metal enrichment, and optionally incorporate galactic 
winds powered by supernovae, as implemented by Springel 
& Hernquist (2003a) . Furthermore, wo follow the growth and 
feedback of BHs based on a new model that combines the 
prescriptions outlined in Springel et al. (2005b) and Sijacki 
& Springel (2006). 

In the following, we briefly summarize the main features 
of the BH model relevant for this study (see Springel et al., 
2005b, for further details), focusing on extensions that per- 
mit us to follow BH growth and feedback in a cosmological 
simulation of structure formation (see also Di Matteo et al., 
2007). We then discuss a new BH feedback prescription at 
low accretion rates, based on a 'bubble' heating scenario 
for representing AGN heating in the radiatively inefHcient 
regime of accretion. 



2.1 Black hole formation and growth 

The BHs in the code are represented by coUisionless sink par- 
ticles, which may accrete gas from their surroundings based 
on a prescribed estimate for the accretion rate. Two BH par- 
ticles are also allowed to merge if they fall within their local 
SPH smoothing lengths and if their relative velocities are 
smaller than the local sound speed. 

In our cosmological simulations of structure formation, 
we assume that low-mass seed BHs are produced sufficiently 
frequently such that any halo above a certain threshold-mass 
contains one such BH at its centre. Whether these seed BHs 
originate in exploding pop-Ill stars, in the collapse of star 
clusters, or are of primordial origin is not important for our 
analysis, but there needs to be a process that produces initial 
seed BHs which can then grow to the masses of supermassive 



BHs by gas accretion in the course of our simulations. For 
definiteness, in most of our simulations we adopt a seed BH 
mass of 10^ h~^MQ and endow all halos with a mass larger 
than 5 x 10^" h~^M(T) with a seed if they do not contain any 
BH already. We identify halos without BHs on the fly during 
a simulation by frequently calling a fast, parallel friends-of- 
friends algorithm that is built into our simulation code. 

State-of-the-art cosmological simulations of structure 
formation reach mass resolutions that are at best of order 
of our initial BH seed mass, while the resolution typically 
reached is still considerably coarser than that. This would 
mean that the initial growth of the BH mass could be signif- 
icantly affected by numerical discreteness effects if the sink 
particle can only swallow full gas particles, as is the case in 
our scheme. In order to avoid that the BH growth and the 
accretion rate estimate is strongly affected by this numeri- 
cal discreteness limitation, we treat the BH mass Mbh as an 
internal degree of freedom of the sink particle. In the begin- 
ning, Mbh may differ from the dynamical mass M^yn of the 
sink particle itself The variable A^bh is integrated smoothly 
in time based on the estimated accretion rate onto the BH, 
while M(jyn increases in discrete steps when the sink particle 
swallows a neighbouring gas particle. The latter process is 
modelled stochastically such that Mdyn tracks Mbh in the 
mean, with small oscillations around it. With this prescrip- 
tion for the BH mass, we can follow the early growth of BHs 
accurately in a sub-resolution fashion even when their mass 
may be smaller than the mass resolution of the simulation, 
while at late times (or in the limit of very good mass res- 
olution) the two mass variables coincide. Of course, in the 
event of a merger of two BH sink particles, both Mbh and 
Mdyn are added together. 

Following Di Matteo et al. (2005), we estimate the ac- 
cretion rate onto a BH particle according to the Bondi- 
Hoyle-Lyttleton formula (Hoyle & Lyttleton, 1939; Bondi 
& Hoyle, 1944; Bondi, 1952) 



Mbh 



4:71 a G'^Miup 



\3/2 



(1) 



where a is a dimensionless parameter, p is the density, Cs 
the sound speed of the gas, and v is the velocity of the BH 
relative to the gas. We account for the possibility that the 
BH accretion has an upper limit given by the Eddington 
rate 

4 TT G Mbh Wp 



MEdd = 



(2) 



er (Tt c 

where rUp is the proton mass, ctt is the Thompson cross- 
section and er is the radiative efflciency, that we assume to 
be 0.1, which is the mean value for the radiatively efficient 
Shakura & Sunyaev (1973) accretion onto a Schwarzschild 
BH. In some of our numerical models we specifically ex- 
plore the imprints of the imposed Eddington limit on the 
BH properties. 

2.2 Black hole feedback 

In the model of Springel et al. (2005b), it is assumed that a 

fixed fraction of the BH bolometric luminosity couples ther- 
mally to the local gas, independent of the accretion rate 
and environment. In our model we extend this BH feedback 
prescription in order to obtain a physically refined model for 
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AGN heating both at high and at low accretion rates. We are 
motivated by the growing theoretical and observational evi- 
dence (Fender et al., 1999; Gallo et al., 2003; Churazov et al., 
2005; Heinz et al., 2005; Croton et al., 2006) that AGN feed- 
back is composed of two modes, analogous to states of X-ray 
binaries. Specifically, at high redshifts and for high accretion 
rates we assume that the bulk of AGN heating is originating 
in the luminous quasar activity. In this regime, BHs accrete 
efficiently and power luminous quasars where only a very 
small fraction of their bolomotric luminosity couples ther- 
mally to the gas. On the other hand, at lower redshifts and 
for BHs accreting at much lower rates than their Edding- 
ton limits, AGN heating proceeds via radiatively inefficient 
feedback in a mostly mechanical form. 

To model the transition between these two accretion 
and feedback modes, we introduce a threshold Xradio = 
Mbh/Msm for the BH accretion rate (BHAR) in Eddington 
units, above which 'quasar heating' is operating, and below 
which we deal with "radio" mode feedback, which we model 
by injecting bubbles into the host galaxy/cluster. Typically, 
we adopt a value of 10~^ for Xradio, and we impose no other 
criterion to distinguish between the two modes of feedback. 

For BHAR that are higher than Xradio, wo parameter- 
ize the feedback as in Springel et al. (2005b), i.e. a small 
fraction of the bolometric luminosity is coupled thermally 
and isotropically to the surrounding gas particles, with an 
amount given by 

Efeed = CfLr = ef MbH C . (3) 

Here ef gives the efficiency of thermal coupling. The value of 
5% adopted here brings the simulated Mbh — o", relation ob- 
tained for remnants of isolated galaxy mergers in agreement 
with observations, as shown by Di Mattco et al. (2005). 

Below Xradio we assume that the accretion periodically 
produces an AGN jet which inflates hot bubbles in the sur- 
rounding gas. We clearly lack the numerical resolution for 
self-consistent ab initio simulations of the detailed physics 
of BH accretion and the involved rclativistic MHD that is 
responsible for the actual jet creation. However, we can nev- 
ertheless try to represent the relevant heating mechanism by 
directly injecting the energy contained in the AGN-inflated 
bubbles into the ICM. To this end we need to link the bub- 
ble properties, like radius, duty cycle and energy content, 
directly with the BH physics. 

Up to now there is no compelling theory that can satis- 
factorily explain AGN jet formation, bubble inflation by the 
jet, and the duty cycle of jet activity. On the other hand, a 
growing body of observational evidence (e.g. Birzan ct al., 
2004; McNamara et al., 2005; Forman et al., 2006; Dunn & 
Fabian, 2006; Fabian et al., 2006) shows that AGN-driven 
bubbles arc present in many systems, at diff'crcnt redshifts 
and over a range of masses, which constrains the duty cycle 
to be of the order of lO'^ — 10* years. However, it is still not 
clear whether and how AGN activity can be influenced by 
properties of the host galaxy like its mass, dynamical state, 
or central gas cooling rate. In light of this uncertainty wc 
propose a simple model for radio feedback where we assume 
that an AGN-driven bubble will be created if a BH has in- 
creased its mass by a certain fraction (5bh = ^Mbh/Mbh. 
Note that with this choice we do not constrain the possible 
duty cycle of the jet itself. We only conjecture that whenever 
the BH increases its mass by 5Mbh, the thermodynamical 
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Table 1. Numerical parameters of our simulations of isolated 
galaxy clusters. The virial masses and radii of the halos, evalu- 
ated at 200 Pcrit, are given in the first two columns. The assumed 

values for the concentration ijararneter are in the third column, 
while the number and the mass of the gas particles is shown in 
the fourth and the fifth column, respectively. The mass of the 
star particles is half that of the gas particles, because we set the 
number of generations of star particles that a gas particle may 
produce to two. Note that there are no parameters for the dark 
matter particles in these run, because we here modelled the dark 
halo with a static NEW potential. Finally, in the last column, the 
gravitational softening length e for the gas and star particles is 
given. 

state of the surrounding gas will be affected by the BH ac- 
tivity, and that this can be represented in the form of the 
bubbles. 

Thus, we relate the energy content of a bubble to the 
BH properties as 

Ebnh = Sm Er C 5MbH , (4) 

where em is the efficiency of mechanical heating by the bub- 
bles. Moreover, we link the bubble radius both to JMbh and 
to the density of the surrounding ICM, in the following way 

„ T3 ( -E^bub/ -Ebub,(l ^ ' 

-Kbub = -Kbub.O -, , \p) 

\PICM/PICM,0 y 

where -Rbub.o, Sbub.o, and picM.o are normalization con- 
stants for the bubble radius, energy content and ambient 
density, respectively. The scaling of the bubble radius is mo- 
tivated by the solutions for the radio cocoon expansion in 
a spherically symmetric case (Scheuer, 1974; Begelman & 
Cioffi, 1989; Heinz et al., 1998). With this parameterization 
for J?bub, we mimic a scenario in which a more powerful 
jet will inffate bigger radio lobes, and where a higher ICM 
density will confine the size of the buoyant bubbles more. Fi- 
nally, the spatial injection of the bubbles is random within 
a sphere with radius twice the bubble radius and centred on 
the BH particle in consideration. 

3 SELF-REGULATED BUBBLE FEEDBACK IN 
ISOLATED HALO SIMULATIONS 

We have first carried out a number of simulations of isolated 
galaxy clusters in order to test our BH feedback model and 
to explore its parameter space. These simulations consist of 
a static NFW dark matter halo (Navarro et al., 1996, 1997) 
with a gas component initially in hydrostatic equilibrium. 
The initial gas density profile has a similar form as the dark 
matter, but with a slightly softened core, as explained in 
more detail in Sijacki & Springel (2006). We construct initial 
conditions for halos with a range of masses (see Table 1 for 
the details of the numerical setup) and evolve them non- 
radiatively for a time 0.25 tnubbie to damp out possible initial 
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Figure 1. Panels (a) and (b) show the time evolution of the BH mass and the accretion rate in Eddington units for a BH that has been 
introduced in the centre of an isolated cluster of mass 10^* /i^^Mq. The different curves and symbols give results for increasing <5bhi 
labelled in the panels. Panel (c) illustrates the SFR as a function of time for the same set of simulations (same colour-coding as panel {d)). 
Finally, panels (d), (e) and (/) give the gas density, mass- weighted temperature and entropy radial profiles for this cluster at t = 0.2tH 
in runs without AGN feedback (blue continuous lines) and with AGN heating, where the different curves are for the increasing value of 
i5bH) a-s indicated on the panels. The vertical dotted lines mark the gravitational softening length and the virial radius, respectively. 
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transients. Then, we "switch on" radiative gas cooling and 
star formation, and also introduce a seed BH particle in the 
halo centre. In the following subsections, we describe the 
results for the subsequent growth of the BHs. 



3.1 Exploring the parameter space 

In this section, we analyse the sensitivity of our model and of 
the resulting BH and host halo properties with respect to the 
parameter choices we have adopted. We begin by considering 
an isolated halo of mass 10^'* ^"^Mq , comparing simulations 
with and without BHs. For these numerical tests, we only 
consider AGN feedback in the form of bubbles, while in all 
full cosmological simulations we will include both modes of 
AGN feedback introduced in Section 2.2. 

First, we vary the threshold for bubble triggering in 
terms of accreted mass, (5bh, from 0.1% to 5%, which will 
affect both the number of bubble events and their mean 
energy. The other two free parameters in the BH feedback 
model, Em and the normalization value for -Rbub, are kept 
fixed in this first series of simulations in order to facilitate 
the comparison. For definiteness, for most of the isolated 
halo simulations we select tm ~ 1, i^bub.o ~ 30/i~^kpc, 
Sbub.o = 5 X lO'^^erg, and picM.o = 10*' /i^Mgkpc"^. Note 
that our choice for these parameters will be slightly different 
in cosmological simulations, as we will discuss in Section 4. 
In most of our numerical experiments we start with a small 
BH seed, typically equal to 10^ /i~^M0, to establish if the 
BH growth is self-regulated and whether it leads to a real- 
istic BH mass when a quasi-stationary state is reached. In 
a number of test simulations we also explore a scenario in 
which a massive BH is already present at the very beginning. 

In panels (a) and (6) of Figure 1 we show the BH mass 
and the accretion rate expressed in Eddington units as a 
function of time, for three different values of 5bh = 0.1%, 
1% and 5%. It can be seen that the BH is initially growing 
rapidly starting from the seed mass of 10^ /i-^M©, and the 
accretion rate onto the BH is quite high, reaching the Ed- 
dington level at f = 0.06 tn. However, after this initial phase 
of rapid growth, AGN feedback starts to reduce the further 
BH growth, because at this point enough heating is sup- 
plied to the surrounding medium, preventing it to quickly 
cool and sink towards the most inner regions. Thus, the 
feedback provided by AGN-driven bubbles reduces the sup- 
ply of gas available for accretion by the central BH. Con- 
sequently, Mbh drops and the mass of the BH shows no 
significant growth with time for a while. However, after a 
certain amount of time which depends on the heating effi- 
ciency, the central cluster gas starts to cool again, causing an 
increase of the BHAR and the triggering of another bubble 
episode. This mechanism establishes a self-regulated cycle, 
in which a balance between the gas cooling rate, the bubble 
heating rate and the residual BHAR is achieved. In panel (b) 
of Figure 1, this cycle of AGN activity can be clearly seen. 
Here the peaks in BHAR correspond to the bubble injection 
events. Note that the higher values of (5bh lead to fewer and 
more powerful AGN outbursts that cause the accretion rate 
to drop to very low values of order of 10~*MEdd- 

The (5bh parameter also affects the properties of the 
ICM. In panels (d), (e) and (/) of Figure 1 we illustrate this 
effect by plotting the radial profiles of gas density, temper- 
ature and entropy, for the different values of iSbh assumed 
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Figure 3. Growth of the BH mass with time (top panel) for 
three clusters of different mass. The middle panel shows how the 
thermal energy contrast, AEbub thermal i depends on time and 
the BH mass in our cluster simulations. The bubble radius as a 
function of time is shown in the bottom panel, for the same set of 
runs. The blue, continuous line and the diamonds are for a cluster 
of mass 10^^ h~^M.Q, the green, dotted line and triangles are for 
a 10^* /i~^Mq halo, while the red, dashed line and the squares 
denote results for a 10"*^^ h~^TsAp, cluster. 
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Figure 2. The panels on the left-hand side show mass-weighted projected temperature maps of three isolated clusters with increasing 
mass: 10^'' h~^M0 (upper-most panel), lO^** h~^MQ (central panel) and 10^^ /i~^Mq (lower-most panel), ft can be seen that the bubbles 
injected in the central cluster regions have mushroom-like morphologies and are uplifting residual cool material from the centre (lower- 
most panel). The middle and right-hand panels give projected pressure and X-ray emission maps, divided by the corresponding map 
smoothed on a scale of 16/i~^kpc. A number of bubble-induced sound waves and weak shocks are clearly visible, and in the case of the 
10^^ fc~^M0 cluster, NW from the cluster centre a mushroom-shaped bubble in the pressure map is clearly visible which corresponds to 
the red blob in the temperature map. 



in the simulations, as indicated in the panels. We also show 
the corresponding profiles for the run without AGN heat- 
ing, denoted by the continuous blue lines. Clearly, there is 
a trend in the simulated profiles with 5bh: more frequent 
and more gentle bubble heating events increase the central 
cluster entropy less than more powerful but rarer feedback 
episodes. A similar result has been found in the simulations 
of Omma & Binney (2004), where a higher jet power was 
heating the cluster for a longer time interval. However, re- 



gardless of the exact value of the 5bh parameter, the cen- 
tral cooling catastrophe is prevented in our model cluster 
in all simulated cases with feedback. Once the feedback be- 
comes effective, the BH mass stops growing when it reaches 
3 — 6 X lO'^ /i~^M0. Note however that the final value of 
the BH mass depends also on the parameter, which we 
have taken to be as high as possible in these test runs. Con- 
sequently, the BH masses obtained here represent a lower 
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Figure 4. Radial profiles of gas density (upper-most panel), 
temperature (central panel) and entropy (lower-most panel) for 
10^^ /i-^M© and lO^'^ h.-'^MQ galaxy clusters at t = 0.2 tn- The 
blue, continuous lines correspond to the simulations without AGN 
heating, while red, dashed lines give the case where AGN feedback 
is present. On every panel, the upper set of profiles are for the 
more massive cluster. Vertical dotted lines mark the gravitational 
softening length and the virial radii of these clusters. 



limit for the actual BH mass that one would expect in a 
cluster of this size. 

In panel (c) of Figure 1, we show how the star forma- 
tion rate (SFR) is affected by AGN-driven bubbles. Without 
AGN heating, the SFR of the lO" ft-^M© isolated cluster 
grows with time to values as high as ~ 25OM0 /yr. This high 
SFR results from the large amounts of cool gas that accu- 
mulate in the central regions, forming a reservoir for intense 
star formation activity. However, in the models with AGN 
feedback, the SFR drops dramatically to values of order of 
~ IMq /yr. For these runs it can be also seen that the SFR 
shows a number of short spikes of somewhat enhanced ac- 
tivity, and these peaks are directly related with episodes of 
increased BHAR which lead to a triggering of bubbles and 
a successive reduction of the SFR. 

We have also tested to which extent the AGN feedback 
is affected by changes in the normalization value for the bub- 
ble radius. As before, we adopt as default values for bubble 
energy and ICM density 5 x lO^^erg and lO*' /i^Mokpc"^, 
respectively, but we vary i?bub,o for different runs. In the 
range of 15/i~^kpc to 30/i~^kpc for iibub.o, our results are 
not altered significantly by modifications of this parame- 
ter. However, for substantially larger values of -Rbub,o, e.g. 
60 /i~^kpc, the heating energy per particle drops to consider- 
ably smaller values, changing the initial rapid growth phase 
of the BH. For large values of -Rbub,o, the AGN feedback is 
less efficient in the early phase of growth such that the BH 
can grow to a somewhat larger mass before a self-regulated 
heating loop is established. We will further discuss the im- 
portance of the bubble radius when we consider halos of 
different mass in Section 3.2. 

Finally, we have numerically explored scenarios where 
an already massive BH is introduced at the very beginning 
of a simulation. For 5bh equal to 0.1%, the BH mass is 
found to increase from its initial value of 5 x IQ^ h'^V^Q 
to 9 X lO^/i~^M0 over a time span equal to one quarter 
of the Hubble time. Most of this BH growth, however, can 
be attributed to the intermittent nature of bubble heating, 
while in the continuous AGN feedback regime, the BH mass 
increases by less then 1% over the same simulated time span. 
Thus, our numerical scheme indeed leads to a stable galaxy 
cluster solution where overcooling in the central parts is pre- 
vented and the BH mass grows to a self-regulated value. 

3.2 AGN heating in halos of different mass 

In this section we examine AGN feedback effects in clus- 
ter halos spanning a range in mass, from 10^^ /i~^Mq to 
10^^ /i~^Mq . We adopt one set of feedback parameters for 
all simulations, namely em, J?bub,o, -Ebub.o, and picM.o, as in 
the previous section, and for 5bh we choose the intermediate 
value of 1%. We are interested in the question whether for 
a fixed set of parameters our model produces satisfactory 
solutions for halos of very different masses. 

Before we delve into a detailed quantitative analysis of 
the cluster simulations, we illustrate the visual morphology 
of simulated AGN-infiated bubbles in Figure 2. For this pur- 
pose, we show projected mass-weighted temperature maps 
(first column), pressure fiuctuation maps (middle column) 
and fiuctuation maps of the X-ray luminosity (last column). 
Different rows are for clusters of increasing mass, from top 
to bottom. We call the maps in the second and the third 
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columns "fluctuation" maps because they have been con- 
structed by dividing the original map by a version that has 
been smoothed on a scale of 16 /i~^kpc. This highlights local 
departures from the mean ICM properties, similar to what is 
known as unsharp masking technique in observational anal- 
ysis. 

A number of interesting features can be seen from the 
panels in Figure 2. First, AGN-driven bubbles have charac- 
teristic mushroom-like and cap-like morphologies. These are 
particularly evident for the 10 ' h'^ Mg halo, where several 
bubbles injected briefly after one another can be noticed. 
Note that the reason why there are so many bubbles in this 
particular map is in part due to the selected time for the 
image. At this moment Mbh happens to grow rapidly, thus 
the bubble duty cycle is rather short. In the last row of Fig- 
ure 2, two bubbles can be seen in the central region of the 
10^^ h-^Ms cluster, one SW from the cluster centre and the 
other NW and more distant. The latter one is uplifting some 
residual cool gas from the cluster centre, which forms a fila- 
mentary structure in the wake of the bubble. Looking at the 
pressure map of this cluster a nicely defined mushroom-like 
bubble can be seen that corresponds to the NW red blob 
in the temperature map. In general, the pressure and X-ray 
luminosity maps show numerous irregularities, sound waves 
and weak shocks that have been generated by the bubbles. 
The ripples in these maps become progressively weaker for 
more massive cluster, which is a consequence of the param- 
eterization adopted in our model, as we will clarify in the 
following paragraphs. 

In the upper- most panel of Figure 3, wc show the mass 
growth of the central BH as a function of time for the dif- 
ferent cluster simulations. It can be seen that halos of in- 
creasing mass exhibit a qualitatively similar behaviour for 
their BH growth: an initial rapid growth phase is followed 
by self-regulated stagnation. However, the starting point of 
the rapid growth is dictated by the gas cooling time, which 
is shortest for the smallest cluster mass considered. There- 
fore, for the 10^'* /i^^Mq halo, the BH starts growing first, 
but it is also the first one to reach stagnation, while the 
BH sitting in the centre of the lO'^^ /i~^Mq cluster waits for 
~ 0.05 t/tn before increasing its mass significantly. At the 
end of the simulated time-span, all three BHs reach equilib- 
rium states and their final masses are: ~ 6.5 x 10® h~^MQ, 
~ 5 X 10'' /i"^Mq, and ~ 1.8 x 10* H'^Mq, respectively. 

In the other two panels of Figure 3, we show how certain 
bubble properties evolve with time for different cluster sim- 
ulations. In the middle panel, AEbub, thermal represents the 
thermal energy contrast of the bubble just before and after 
the energy was injected into it. AEbub, thermal is reflecting 
mainly the mass growth of the BH, since it determines how 
much energy will be thermally coupled to the bubbles. It 
can be seen that AEbub, thermal initially grows significantly 
with time, but when the BH growth saturates, it scatters 
around a constant value. Also, there is a systematic trend of 
AEbub, thermal with cluster mass, with smaller mass clusters 
showing a larger energy contrast for the bubbles. This is pri- 
marily due to the smaller radius of bubbles in less massive 
clusters, as can be seen in the bottom panel of Figure 3. 
Given that we have adopted the same set of parameters 
for our different cluster simulations, it follows from Eqn. (5) 
that BHs of smaller mass will generate smaller bubbles. This 
explains why Rhuh is growing in time and why it is bigger 



for more massive halos that harbour larger BHs. This conse- 
quence of our model has also repercussions for the presence 
and intensity of sound waves generated by the inflation of 
bubbles, as seen in the maps of Figure 2. Smaller bubbles 
with a higher energy contrast relative to the surrounding 
ICM are more efficient in generating weak shocks and sound 
waves than large bubbles with low AEbub, thermal. 

In Figure 4, we show radial profiles of gas density, 
temperature and entropy for our isolated clusters of mass 
lO"ft~iM0 and 10^'"^ ft-^M©. The radial profiles of the 
10" h-^MQ cluster can be found in Figure 1 with the same 
colour-coding. The blue, continuous lines denote the runs 
without AGN, while the red, dashed curves are for the sim- 
ulations with AGN heating. The upper set of profiles in each 
panel arc for the more massive cluster. The results show that 
for a fixed choice of parameters in our AGN heating model 
the ICM properties of clusters over a range of masses appeax 
satisfactory: (i) the gas density is somewhat suppressed in 
central regions; (ii) very cool gas in the innermost regions 
of clusters is absent, and at the same time the AGN heating 
is not increasing the ICM temperature too much; (iii) the 
gas entropy is boosted in central regions, but without gen- 
erating entropy inversions that are typically not observed in 
real systems. Together with realistic BH masses and BHAR, 
these promising features of the simulated ICM properties 
encourage us to proceed studying our model using fully cos- 
mological simulations of galaxy cluster formation, which we 
consider in the next section. 



4 COSMOLOGICAL SIMULATIONS OF AGN 
FEEDBACK IN GALAXY CLUSTERS 

We now turn to self-consistent cosmological simulations of 
galaxy cluster formation, with the aim to understand how 
AGN feedback infiucnces these objects at different cosmolog- 
ical epochs and how this influence depends on the mass and 
the evolutionary history of clusters. We focus our analysis 
on two galaxy cluster simulations that have quite different 
merging histories and different present-day masses. The clus- 
ters have been selected from a cosmological ACDM model 
with a box size of 479fe"^Mpc (Yoshida et al, 2001; Jenkins 
et al., 2001), and were prepared by Dolag (2004) for resimu- 
lation at higher resolution using the Zoomed Initial Condi- 
tion (ZIC) technique (Tormen et al., 1997). The cosmological 
parameters of the simulations are those of a ACDM concor- 
dance cosmology with f2m = 0.3, Oa = 0.7, Ob = 0.04, 
as = 0.9, and Hubble constant Ho = 70kms~^Mpc~^ at 
the present epoch. In respect to the isolated halos sim- 
ulations we have adopted somewhat different parameters 
of the BH feedback model, with the normalization values 
for i;bub,o = lO'^^erg, and picM,o = lO'' M© kpc"^ being 
lower, with Jbh = 0.01%, and with a more realistic me- 
chanical feedback efficiency of tm = 0.20, while the other 
parameters were kept exactly the same. 

The primary numerical parameters of our sinmlations 
are summarized in Table 2, and the main physical properties 
at = of the formed clusters are listed in Table 3. For both 
simulated galaxy clusters, the smaller one ("g676") and the 
more massive one ( "gl" ) , the virial radius, virial mass, total 
gas mass and stellar mass have been measured in runs with- 
out AGN feedback (labelled "csf") and when it is included 
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Figure 5. Projected gas density map of tlie g676 galaxy cluster 
simulation at 2 = 0, subject to AGN heating. Black dots mark 
the positions of BH particles that are more massive than 1.5 X 
10^ /i~^Mq. It can be seen that the BHs are residing not only in 
the centre of the most massive galaxy cluster, but also in almost 
all smaller halos that are visible on the map. Interestingly, there 
is a concentration of BHs in a filament that is extending from 
the central cluster towards NW, due to the ongoing structure 
formation there. 

(labelled "csfbh"). Additionally, the mass- and emission- 
weighted temperature. X-ray luminosity and total SFR are 
given in the 6th to the 9th column of Table 3, respectively. 
Looking at the gas properties listed in Table 3 it can be seen 
that AGN feedback does not significantly alter the global 
gas temperature, or the X-ray luminosity. However, there 
are still important local changes of these properties when 
AGN heating is included, as we will discuss in more detail 
in Section 4.3. On the other hand, the global stellar proper- 
ties are strongly affected by AGN feedback: the total stellar 
mass is reduced considerably, as well as the total SFR, while 
the total gas mass is increased. This implies that the rela- 
tive amount of gas versus stars is a function of AGN heating 
efficiency. 

4.1 Black hole growth in clusters 

In Figure 5, we show a projected gas density map of the 
g676 galaxy cluster simulation at z = 0, with the positions 
of BHs more massive than 1.5 x 10^ /i~^M0 overlaid as black 
dots. The density map is 20/i~^Mpc on a side and is cen- 
tred on the most massive cluster galaxy. Interestingly, the 
surrounding large-scale structure shows many smaller halos 
with embedded BHs, many of them residing in a filamentary 
region NW of the central object. In fact, most halos are har- 
bouring a central, massive BH. However, the most massive 
BH among them is sitting at the centre of the biggest halo, 
and this is also true at higher redshifts. 

The mass as a function of time of this most massive 



BH is shown in Figure 6. To produce this measurement, we 
have first constructed the full merger tree of all BHs formed 
in the simulated volume, allowing us to extract the main 
progenitor trunk corresponding to the most massive BH at 
2 = 0, which is plotted as a thick continuous blue line. Also, 
we show the mass growth of the secondary progenitors that 
at the moment of their last merger have a mass greater than 
5 X 10^ h~^yiQ. It can be seen that the most massive BH at 
2 = is not the first one to form, but it merges with two BHs 
that formed somewhat earlier and that are more massive for 
z > 3.5. This is similar to our findings in Di Matteo et al. 
(2007), where we showed that the most massive black holes 
found at high redshift in a cosmological volume need not 
necessarily be the most massive ones at low redshift. The 
thick, dashed blue line in Figure 6 represents the cumulative 
mass of all secondary progenitors of the most massive BH 
at 2 = 0. This shows that mergers contribute up to 45% to 
the final mass of the BH that defines the trunk, while the 
remaining 55% are due to gas accretion by the BH on the 
trunk itself. Note however that ultimately all the BH mass is 
built-up by gas accretion (the mass of the initial seed BHs is 
negligible), but the part we here attributed to mergers was 
built up by gas accretion in other systems. 

We also note that the mass growth via mergers is more 
important at higher redshifts, with the last significant merg- 
ing event occurring at z ~ 0.7. The merger events of BHs 
are mainly driven by the merging history of the host galaxy 
clusters. Indeed, there is a quite close correspondence be- 
tween the evolution of the main progenitor of the most mas- 
sive BH and the main progenitor of the biggest halo in the 
simulation. In particular, the main progenitor of the g676 
galaxy cluster undergoes a number of mergers between red- 
shifts 3.5 < z < 2.8, when also the BH grows significantly 
via mergers with other BHs. Then at z ~ 0.75, the g676 
cluster experiences its last major merger, visible as a jump 
in the main trunk in Figure 6. At lower redshifts, the cluster 
becomes fairly relaxed and isolated, such that its central BH 
grows mostly by accreting gas that is cooling off from the 
hot cluster atmosphere. As we will discuss in Section 6 this 
mode of accretion (at low accretion rates) does however not 
contribute significantly to the overall BH mass growth. 

Now we consider in more detail the mass growth and 
the accretion rate, not only of the most massive BH, but of 
all the BHs belonging to the simulated volume. Due to the 
better statistics provided by the gl galaxy cluster simula- 
tion, we focus on this simulation, noting that the results for 
the g676 cluster are very similar. In Figure 7, we plot the 
BHAR, expressed in Eddington units, as a function of the 
BH mass. The colour-coding denotes the redshift, over the 
redshift interval z = 4.3 to z = considered here. A number 
of interesting features can be noticed in this plot: (i) at a 
given BH mass, the accretion rate is decreasing with redshift; 
(ii) for intermediate redshifts, 2 < 2 < 3, the BHAR is low- 
est for intermediate mass BHs; (iii) at still lower redshifts, 
the BHAR for intermediate masses keeps falling strongly, 
while it stays comparatively high for low-mass BHs as well 
as for very massive BHs (iv) at the present epoch there is 
a large population of low-mass BHs that is still accreting 
efficiently. These features are qualitatively consistent with 
observational findings both from optical and X-rays surveys 
(e.g. Steffen et al., 2003; Ueda et al., 2003; Heckman et al., 
2004; Barger et al., 2005; Hasinger et al., 2005) that indi- 
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Simulation 






muM [h ^Mq] 


rrtgas [h ] 






e [/i-^kpc] 


g676 


314518 


314518 


1.13 X 10^ 


0.17 X 10^ 


60 





5.0 


gl 


4937886 


4937886 


1.13 X 10^ 


0.17 X 10^ 


60 





5.0 



Table 2. Numerical parameters of the cosmological galaxy cluster simulations analysed in this study. The values listed from the second 
to the fifth column refer to the number and to the mass of high resolution dark matter particles and of gas particles. Note that the actual 
values of Afgaa and mgas vary in time due to star formation. The last three columns give the initial and final redshifts of the runs, and 
the gravitational softening length e. 
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R200 


M200 


A^gas,200 


Afatars,200 


7mw 


Tew 
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[h-lkpc] 
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[/i-iM©] 


[K] 


[K] 


[ergs-i ] 




g676_csf 


1176 


1.13 X 10^'' 


9.3 X 10^2 


4.7 X 10^2 


1.4 X 10^ 


2.6 X 10^ 


1.6 X 10*^ 


51 


g676_csfbh 


1165 


1.10 X 10^'' 


1.1 X 10^^ 


1.4 X 10^2 


1.3 X 10^ 


2.4 X 10^ 


1.0 X 10**^ 


1 


gl_csf 


2857 


1.63 X 10^5 


1.4 X 10" 


6.3 X 10^^ 


7.3 X 10^ 


1.3 X 10** 


1.0 X 10-*5 


742 


gl_csfbh 


2859 


1.63 X 10^5 


1.7 X lO^-* 


2.5 X 10" 


7.3 X 10^ 


1.3 X 10** 


1.0 X 10^5 


144 



Table 3. Physical properties of our sample of simulated galaxy clusters at ^ = 0, selected with a mean overdensity of 200pc. For two 
different galaxy clusters, as labelled in the first column, the cluster radius, total mass, gas mass and stellar mass are given. Also, the 
mass— and emission-weighted gas temperatures, X-ray luminosity and total SFR are listed, in the 6th to 9th column, respectively. Note 
that the values give results both for the simulations with cooling and star formation only (denoted with "csf"), and for runs performed 
with AGN heating (labelled "csfbh"). 
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Figure 6. Merger tree of the most massive BH at z = in the 
g676 galaxy cluster simulation. The thick blue line is showing 
the mass growth of the BH's first progenitor as a function of 
redshift, while the thin lines of different colour represent the evo- 
lutionary history of secondary progenitors, until the merger with 
the main progenitor occurs. For clarity, only secondary progeni- 
tors that at the moment of the last merger have a mass greater 
than 5 X 10^/i~^Mq are shown. Vertical dotted lines indicate 
when merger events between the BHs in consideration occur. The 
dashed, thick blue line represents the cumulative mass of all sec- 
ondary progenitors. 



cate that BHs are growing in a so-called 'anti-hierarchical' 
manner. In Section 5, we will come back to this issue, and 
discuss how the BHAR depends on the BH mass in a galactic 
environment. 

At low redshifts, the increase of the BHAR with mass 
at the high BH mass end is driven by BHs that sit at the 
centres of massive galaxy clusters. These BHs are fed by 
gas from the hot cluster atmosphere, which would develop 
a strong cooling flow without the periodic heating by the 
AGN feedback. Interestingly, we also find that our BHAR 
of massive BHs ai z — agrees very well with a recent 
estimate of the Bondi accretion rate for a sample of X-ray 
luminous elliptical galaxies by Allen et al. (2006). 

However, we need to emphasize that the simulations 
performed in this study were not designed to address the 
problem of BH formation and growth in the very early uni- 
verse, for which they lack the required resolution and vol- 
ume. This becomes also apparent in the plot of Figure 7, 
where the upper right corner is not populated. This means 
that these simulations cannot account for the existence of 
supermassive BHs of mass ~ 10^ Mq already at z = 6, as 
inferred from the luminosity of high redshift quasars (Fan 
et al., 2001). Using high-resolution compound galaxy models 
to populate a merger tree measured from a zoom simulation 
of a protocluster region, Li et al. (2006) have however re- 
cently demonstrated that at least in principle such early su- 
permassive BHs can grow from accretion in gas-rich mergers 
at high redshift. 

4.2 Heating the cluster outskirts 

In Figure 8, we examine the radial dependence of AGN heat- 
ing in our simulated clusters for different redshift bins. To 
this end, we evaluate the AGN luminosity at the given epoch 
and plot it as a function of distance from the centre of the 
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Figure 7. BHAR in Eddington units as a function of BH mass, for 
all BHs belonging to the gl galaxy cluster simulation. The colour- 
coding expresses the redshift, and the redshift interval considered 
here ranges from 4.3 to 0.0. For a given BH mass, the accretion 
rate is highest at early cosmic times. At lower redshifts, there is 
an upturn in the BHAR at the massive BH end, corresponding 
to BHs that are fed by gas that cools off from the atmosphere of 
the cluster and its massive progenitor systems. 



most massive cluster in the simulation. Depending on the 
BHAR, we decide whether the BH is in the "quasar" or "ra- 
dio" mode, which we denote with star symbols or circles, 
respectively (see upper panel). In the lower panel, we plot 
the total AGN luminosity outside of a given radius, regard- 
less in which mode the BHs accrete. It can be seen that the 
AGN heating at all redshifts considered and regardless of 
the feedback mode is most important in the cluster centre. 
However, at early times, and in particular for 1.5 < z < 4.3, 
BHs that happen to reside in cluster outskirts during this 
time could provide an important additional source of ICM 
heating. These BHs are essentially all in the quasar phase 
where they accrete gas efficiently. Most of these BHs reside 
in satellite halos that are entering the most massive system 
for the first time, and thus possibly are preheating the gas 
of smaller halos prior to and during the merger with the 
massive cluster. These findings and the spatial distribution 
of AGN in galaxy clusters, illustrated in Figure 5, are in 
a good agreement with observational evidence (e.g. Cappi 
et al., 2001; Cappelluti et al., 2005; Ruderman & Ebeling, 
2005; Martini et al., 2006) for the presence of AGN in galaxy 
cluster environments. At low redshifts, though, the AGN sit- 
ting at the cluster centre appears to be the dominant source 
of heating, as we confirm in Section 4.3. Interestingly, the 
feedback luminosity of the central AGN, regardless of its ac- 
cretion regime, peaks at the similar value for the whole time 
span considered. 
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Figure 8. In the upper panel, the mechanical and quasar AGN 
luminosities are plotted as a function of distance from the most 
massive halo in the g676 galaxy cluster simulation. Depending on 
the BHAR, a BH is assumed to be either in a "quasar" or in a 
"radio" mode, and the BH luminosity is calculated accordingly 
and denoted with different symbols. The colour-coding is accord- 
ing to redshift bins, as indicated in the lower panel. The total 
AGN luminosity outside a given radius, regardless of the feed- 
back mode, is shown in the lower panel. It can be seen that the 
highest AGN luminosities always correspond to BHs sitting in the 
centre of the most massive clusters, and that at higher redshifts 
(especially 1.5 < 2 < 4.3) heating from the quasars in cluster 
outskirts can be important. 
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Figure 9. Radial profiles of gas density (upper row), mass-weighted temperature (central row) and entropy (lower row) of the g676 
galaxy cluster at 2 = 0.5, 0.2 and 0.0, respectively. Continuous blue lines illustrate the case without AGN heating, while dashed red lines 
are for simulations where AGN feedback is included. The run where the Eddington limit was not imposed on the BHAR is shown with a 
dotted green line. The vertical dotted lines denote the gravitational softening and the virial radius of this cluster. The dash-dotted lines 
in the central row of panels show the slope of the central temperature profiles of the cool core clusters found by Sanderson et al. (2006), 
while the dash-dotted lines in the lower row illustrate the entropy scaling with radius, i.e. S oc r^'^. 
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Figure 10. Radial profiles of mass-weighted gas mctallicity and stellar metallicity of the g676 (left-hand panel) and gl (right-hand 
panel) galaxy cluster simulations at z = 0. Continuous lines show the metallicity profiles without AGN heating, while dashed lines are 
for the runs with AGN feedback. The arrows illustrate how the ICM metallicity is affected by the BHs. The stellar metallicity is reduced 
at all radii, while the gas metallicity shows a tilt in the profile, decreasing in the inner regions and increasing for r > 30/i^^kpc, while 
still exhibiting a residual gradient. 



4.3 The impact of AGN heating on the ICM 

We here analyse how AGN feedback affects the intracluster 
medium and the stellar properties of the cluster galaxies. In 
Figure 9, we show radial profiles of gas density, temperature 
and entropy of the g676 cluster at three different epochs, 
both without (continuous blue lines) and with AGN heating 
(dashed red lines). Additionally, we have explored a model 
in which the BHAR was not limited to the Eddington rate, 
which is drawn with a dotted green line. It can be seen that 
at all redshifts considered, the central gas density is reduced 
significantly as a result of AGN feedback. We also observe 
that the gas temperature in the central regions is fluctuating 
around the "isothermal" value, sometimes for short time in- 
tervals increasing towards the centre, but most of the time 
being roughly constant or decreasing. These trends in gas 
temperature are due to the periodic nature of the bubble 
feedback. If a very energetic bubble is injected in the in- 
nermost cluster regions, the gas temperature increases for a 
short time, and reduces the BHAR considerably. After some 
time has elapsed, the gas will begin to cool again such that 
the central gas temperature will gradually start to decrease 
towards the centre. When enough gas has cooled off from 
the hot cluster atmosphere and has become available for 
accretion onto the supermassive BH, a new bubble will be 
triggered, establishing a self-regulated loop for the growth 
of the BH and the heating of the ICM. In some sense, this 
feedback loop acts as a thermometer for the ICM, preventing 
it to develop strong cooling flows. 

On top of our simulated gas temperature profiles, we 
plot the slopes of the central temperature profiles of the 
cool core clusters recently found by Sanderson et al. (2006) , 
showing a very encouraging agreement. Indeed, a number 



of observational studies (e.g. De Grandi & Molendi, 2002; 
Vikhlinin et al., 2005; Dunn & Fabian, 2006; Pratt et al., 
2007) have shown that the central temperature profiles of 
relaxed, cool core clusters decrease with decreasing radius, 
while at the same time these systems require an additional 
source of heating in order to avoid excessive cooling fiows. 
Theoretically, it is not trivial to explain this observed fea- 
ture of the central cluster temperature - a self-regulated 
feedback mechanism is apparently needed which is sensitive 
to the local properties of the intracluster gas. Our model 
contains such a mechanism and represents to our knowledge 
the first successful attempt to simulate such a process in a 
cosmological environment. 

As a result of the non-gravitational bubble heating, the 
ICM entropy is increased in the central regions, but main- 
tains a monotonically increasing radial profile. Given the 
temperature of this cluster, the overall shape and the nor- 
malization of the simulated entropy profile are consistent 
with observational findings (e.g. Pratt & Arnaud, 2003; Pon- 
man et al., 2003). Observationally, it is found that the en- 
tropy scales with radius roughly as 5 oc r^'^ in the cluster 
outskirts, while departures from this scaling are seen in the 
central regions, out to r ~ 0.1 — 0.2 f?200. As can be noticed 
from the bottom row of Figure 9, our simulated entropy pro- 
files depart from this scaling out to somewhat larger radii, 
and this is also found for the 'gl' cluster simulation which 
has quite similar radial profiles. However, to decide whether 
this constitutes a real discrepancy between the simulations 
and observations, a larger set of simulated galaxy clusters is 
required, which we plan to compile in forthcoming work. 

Looking at the radial gas profiles for r > 400/i~^kpc it 
is clear that the AGN heating does not significantly change 
the ICM properties at larger radii. Thus, heating of the clus- 
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Figure 11. Stellar density profile of the g676 (left-hand panel) and gl (right-hand panel) galaxy cluster simulations at 2 = 0. Continuous 
blue lines show the case without AGN heating, while dashed red lines are for simulations where AGN feedback is included. It can be 
seen that the stellar density is decreased at all radii, and in the case of the gl cluster, the central stellar density develops a core when 
AGN feedback is operating. 



ter outskirts does not appear to be relevant for z < 0.5, 
confirming our findings in Section 4.2. Moreover, the radial 
profiles of our simulated galaxy cluster at z = 0.5 do not 
differ significantly from the ones at z = 0. This is consistent 
with observational findings (Bauer et al., 2005), where al- 
ready formed cool core clusters are seen at 2; ~ 0.15 — 0.4, 
with similar properties to the local population. 

In Figure 10, we show radial profiles of gas and stel- 
lar metallicity, in Solar units, for the g676 and gl cluster 
simulations at 2 = 0. The continuous lines show the metal- 
licity profiles without AGN feedback, while the dashed lines 
illustrate how these profiles change when our BH model is 
"switched-on" . The stellar metallicity is reduced at all radii, 
while the gas metallicity shows a tilt: it is reduced in the 
innermost regions, for r < 30/i~^kpc, while it is increased 
for larger radii. These changes in the metallicity gradients 
are caused by two effects: first, due to the AGN feedback 
the total number of stars is reduced, not only in the cen- 
tral regions, but all over the cluster, as illustrated in Fig- 
ure 11; second, some of the metals accumulated in dense, 
star-forming regions are expelled into the hot ICM. This 
mechanism drives the tilt in the gas metallicity, increasing 
the metal content of the hot intracluster gas. 

The effect of AGN heating on the metal production and 
mixing in clusters brings the simulated gradients into a much 
better agreement with observational results. However, the 
simulated gas metallicity in the cluster outskirts appears still 
too low in comparison with the metallicity levels of cool core 
clusters at similar radii (e.g. De Grandi & Molendi, 2001). 
This discrepancy may in part be attributed to the fact that 
in these simulations we have not included feedback effects 
from galactic winds and outflows powered by star forma- 
tion. As we have shown in Sijacki & Springel (2006), galactic 



winds help spreading metals throughout the cluster environ- 
ment, especially into the outer parts. Another restriction of 
our model for metal enrichment lies in its highly simplified 
treatment of supernovae, where we neglect the time delay of 
supernovae type la. More sophisticated models of chemical 
enrichment (e.g. Tornatore et al., 2004; Scannapieco et al., 
2005; Kobayashi et al., 2007) that take this into account 
might therefore help as well. 

We now turn to the effects of AGN feedback on the 
stellar properties of the simulated clusters. In Figure 11, it 
can be seen that the stellar density is reduced at all radii in 
both clusters, as we have anticipated above. This reduction 
is caused not only by the activity of the central cluster BH, 
but also by other BHs that are harboured at the centres of 
individual galaxies throughout the cosmic evolution. As we 
will discuss in more detail in Section 5.2, when a BH at the 
centre of a given galaxy becomes massive enough to influence 
its host, it acts on the stellar properties by reducing both 
the instantaneous SFR as well as the integrated total stellar 
mass by a significant amount. 

In Table 4, we list a number of properties we have mea- 
sured for the central cluster galaxy. For the identification 
of the cD galaxy, we have simply taken its radius to be 
20/i~^kpc, based on the apparent break in the stellar den- 
sity profile at this radius. A number of interesting results 
can be noticed from the values in the table: due to the AGN 
feedback, the stellar mass of the cD galaxy is significantly 
reduced, as well as the central total gas mass. There is some 
residual cold gas in the cD galaxy, but star formation is com- 
pletely quenched. This has an immediate and important con- 
sequence for the photometric colours of the cD galaxies. In 
order to estimate the colours, we have used the stellar pop- 
ulation synthesis models of Bruzual & Chariot (2003) and 
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Table 4. Properties of the cD galaxy of the g676 and gl galaxy clusters at z = 0. The listed values refer both to simulations with cooling 
and star formation only (denoted with "csf"), and to the runs where AGN heating was included (labelled "csfbh"). From the second to 
the fourth column, we give the stellar mass, the gas mass and the cold gas mass with T < 1 keV of the cD galaxy. The star formation 
rate of the cD galaxy is found in the fifth column, while the Mr magnitude and the u — r and g — r colours are given in the last three 
columns, respectively. 
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Table 5. Numerical parameters of our cosmological simulations in periodic boxes. The values listed in the second to the fifth column 
refer to the size of the box, to the number of gas and dark matter particles and to their masses. The last three columns give the initial 
and final redshifts of the runs, and their gravitational softening lengths e. 



computed the rest-frame magnitudes in the SDSS bands, 
assuming Solar metaUicity and a Chabrier initial mass frmc- 
tion. Both cD galaxies in our cluster simulations become 
much less luminous in the r, u and p-bands as a result of 
AGN feedback, such that the u — r and g — r colour indices 
are substantially increased. This reflects the fact that the 
stellar populations of the cD galaxies have become much 
older in our BH model, with very little recent star forma- 
tion. With these properties, the simulated cD gala^cies are in 
quite good agreement with observational findings (e.g. von 
dcr Linden ct al., 2007), in terms of their stellar mass, their 
SFR and their colours. 

Finally, we briefly discuss whether the Eddington limit 
that we have imposed on the BHAR has a significant ef- 
fect on the BH growth and the feedback in clusters. For 
this purpose we have rerun the g676 galaxy cluster simular 
tion with exactly the same feedback parameters, but relax- 
ing the maximum accretion rate assumption. In Figure 9, 
we have included results for the radial profiles of gas den- 
sity, temperature and entropy in this run as dotted green 
lines. It can be seen that the ICM properties are modified 
slightly more strongly in this case with respect to our default 
model, i.e. the AGN feedback effects appear to be somewhat 
stronger. A closer look at the detailed mass growth and ac- 
cretion rate history of the BH as function of redshift shows 
that the feedback is regulating the BHAR even without an 
upper limit in the form of the Eddington rate. In fact, the 
total mass of all BHs at 2 = in the simulated volume 
increases only by roughly 20% if the Eddington limit is dis- 
regarded. The most mgissive BH in the cluster centre shows 
an increase of its mass of the same order. This additional 
growth of the BH mass mostly comes from brief episodes of 
very high accretion rate (typically reaching a few times up to 
ten times MBdd), which are quickly terminated by the onset 
of AGN heating. Without an explicit Eddington limit, a bit 



more growth can happen before the 'suicidal' AGN activity 
shuts off the accretion, but the limit itself is not required 
for establishing a stable self-regulation loop. The Eddington 
limit therefore does not seem to pose an important restric- 
tion on the growth of BHs in clusters. However, this may 
still be quite different in the early phases of BH accretion 
at high redshifts, when the BH mass is so low that feedback 
effects are weak. Then the Eddington limit sets the short- 
est growth timescale that can be realized, which makes it a 
challenge to grow BHs quickly enough to the observation- 
ally inferred large masses already at very high redshifts. We 
plan to address this interesting problem in full cosmological 
simulations in future work. 



5 SIMULATIONS OF GALAXY FORMATION 
WITH AGN FEEDBACK 

In this section we analyse BH growth and feedback in sim- 
ulations of cosmic structure formation in homogeneously 
sampled volumes. We are mainly interested in the question 
whether our numerical model for a two-mode AGN growth 
produces realistic results for a range of object meisses, from 
the scales of isolated galaxies to the ones of massive galaxy 
clusters. In Section 4, we have already confirmed that the 
BH model works well for clusters of galaxies, where it in fact 
drastically improves the properties of simulated galaxy clus- 
ters with respect to observations, and yields at the same time 
plausible evolutionary histories for the BH masses and the 
accretion rates. Here, we want to see whether this success 
also carries over to the scale of individual galaxies, where 
the "quasar" mode of feedback will be dominant. In par- 
ticular, we want to establish whether our unified model for 
AGN feedback maintains and extends the successes we have 
found in Di Matteo et al. (2007) for the high-redshift galaxy 



© 0000 RAS, MNRAS 000, 000-000 



A unified model for AGN feedback in cosmological simulations of structure formation 17 



u I ' ' ' ' I 



o 

Q. 



X 
Cl 



"1 — I — I — I — I — I — I — I — I — I — I — I — I — I I I 



L,,, = 25 h 'Mpc 




o 

Q. 



10-^ r 



© 



o 

X 



10"=" r 



10' 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

L,„, = 25 h'Mpc : 




2x176' 
2 X 256' 
2 X 384' 



1 2 3 4 5 6 7 8 

z 



Figure 12. Rcdshift evolution of the comoving BH mass density in a cosmological box of 25/i~^Mpc on a side (left-hand panel). The 
three different curves show runs with increasing resolution as labelled in the panel. The right-hand panel gives the BHAR (continuous 
lines) and SFR (dashed lines) densities as a function of redshift for the same set of runs. For plotting purposes, psFR has been rescaled 
by a factor of 10'^ to fit onto the same scale. It can be seen that numerical convergence is achieved at low redshift, where BHAR and 
SFR densities of the different runs asymptotically approach each other with increasing resolution. 



population in a model that only accounts for the "quasar" 
mode. 

We have performed hydrodynamical simulations of a 
cosmological box with 25 h~^Mpc on a side, at three differ- 
ent mass resolutions, as summarized in Table 5. The small 
box-size allows us to study comparatively small galaxies 
with stellar masses in the range of ~ 1O*M0 to ~ lO^^M©. 
However, since the fundamental mode of this small box be- 
comes non-linear at around z = 1, we had to stop the sim- 
ulations at this redshift, since they would not be represen- 
tative any more for the low-redshift universe. For the initial 
conditions, we have used a 'glass' as unperturbed particle 
load and the power spectrum fit of Eisenstein & Hu (1999) 
for imposing initial perturbations with WMAP-3 cosmolog- 
ical parameters (Spergel et al., 2007). In particular, we have 
adopted a ACDM cosmology with = 0.26, fib = 0.044, 
r^A = 0.74, (jg = 0.75, and Ho = 71kms"^Mpc~^ at the 
present epoch, and a primordial spectral index of = 0.938. 

For all three mass resolutions R1-R3 that we considered, 
we have computed two runs, one only with cooling and star 
formation and the other with the BH model included as well. 
The parameters of the BH model were exactly the same as 
the ones adopted in our cosmological simulations of galaxy 
cluster formation. For our intermediate resolution case R2, 
we have carried out additional simulations to study the in- 
fluence of further parameters. This includes a run where we 
used the same cosmological parameters as in Section 4 in 
order to gauge the importance of the cosmological model. 
We have also explored the relative importance of the bub- 
ble feedback on galactic scales relative to the quasar feed- 
back, that we will discuss more in detail in Section 6. Fi- 
nally, we have studied the influence of including a model for 
supernova-driven galactic winds with velocity v ~ 480km/s, 



which can be especially important in low mass systems, as 
we will discuss below. 



5.1 Black hole growth 

In Figure 12, we consider the comoving BH mass density 
(left-hand panel) and BHAR density (right-hand panel) as 
a function of cosmic time. For comparison, we also show 
the redshift evolution of the SFR density (right-hand panel, 
dashed lines), which we have rescaled by a constant fac- 
tor of 10^ to put the curves onto the same plot. It can 
be seen that the BHAR and SFR densities approach each 
other asymptotically towards low redshifts when the numer- 
ical resolution is increased. The fact that numerical conver- 
gence can be more easily achieved at low than at high red- 
shift is not surprising in light of the hierarchical growth of 
structure. At high redshift, most of the small mass systems 
that undergo star formation and BH accretion in our high- 
est resolution simulation are unresolved in our low resolution 
simulation, which hence underpredicts the BHAR and SFR 
density. However, towards later times, the mass scale where 
most of the accretion happens shifts to larger mass objects, 
which can be resolved in all of our simulations, such that ap- 
proximate convergence can be reached even with moderate 
resolution. This behaviour is also similar to the one seen by 
Springel & Hernquist (2003b) in a comprehensive simulation 
study of the cosmic star formation rate history. We consider 
it highly encouraging that we find such a behaviour for our 
BH growth model as well, which shows that the model is 
numerically well posed and can be meaningfully applied to 
cosmological volumes with their comparatively poor resolu- 
tion per galaxy. 
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Figure 13. Left-hand panel: BHAR in Eddington units as a function of BH mass, for all BHs belonging to the R2 cosmological simulation. 
The colour-coding expresses the redshift, and the redshift interval considered is from 5.0 to 1.0. Right-hand panel: Characteristic growth 
time of the BHs as a function of their mass, colour-coded according to redshift. The growth time has been defined as the ratio of the BH 
mass and the BHAR at the given epoch, normalized to the Hubble time at z = 0. 



Nevertheless, it is clear that the numerical resolution 
in our cosmological runs remains too poor to correctly re- 
solve the internal structure of galaxies, in particular the thin 
gaseous and stellar disks in spiral galaxies. This presum- 
ably afltects the detailed time evolution of star formation and 
black hole accretion in galaxy mergers. For example, our cos- 
mological runs will largely blend over the time delay we see 
in high-resolution galaxy mergers between BH activity and 
the peak of the star burst. But the final BH mass reached 
in a merger is comparatively insensitive to these details and 
probes primarily the depth of the gravitational potential of 
the forming spheroid, which our simulations can still repre- 
sent quite well. While numerical resolution is clearly an im- 
portant Umitation in our cosmological results, we therefore 
believe they are reasonably accurate for the basic properties 
of the BH population. 

Our estimated BH mass density at 2 = 1 in our highest 
resolution run is 2.5 x lO^M© Mpc""^. This is in very good 
agreement with a number of observational estimates, based 
both on optical and X-ray data, which typically find val- 
ues in the range of 2 - 6 X IO^Mq Mpc"^ (e.g. Fabian & 
Iwasawa, 1999; Merritt & Ferrarese, 2001; Yu & Tremaine, 
2002; Cowie et al., 2003; Graham et al., 2007). The SFR 
density shows a peak at z ~ 3, while the BHAR density 
peaks at somewhat lower redshift, i.e. z ~ 2.5. The BHAR 
and SFR then decline to lower redshifts, with the SFR den- 
sity decreasing slightly more steeply. If we extrapolate the 
BHAR and SFR densities to z = 0, the ratio of the two is of 
order of 2 x 10'^, a value that is roughly in agreement with 
estimates from the local population of galaxies (Heckman 
et al., 2004). In contrast, at high redshift, the BHAR den- 
sity is rising more steeply than the SFR density, in broad 
agreement with the estimates by Merloni et al. (2004). 

However, Hopkins et al. (2007b) find that the quasar 
luminosity density declines more steeply for z > 2 than in 



our simulations. It is quite plausible that this discrepancy 
at high redshift is due to the limited numerical resolution 
available in our cosmological simulations, and in particu- 
lar, their inability to correctly reproduce galactic disks of 
the right sizes, as discussed above. While the bulk of the 
BH accretion happens during major merger events in our 
simulations, we expect that there is a significant contribu- 
tion from minor merger episodes or secular processes, which 
boost the BHAR at high z and induce a less steep rise. If 
the galaxies could be fully spatially resolved, the gas would 
instead not be able to easily reach the galactic centres in 
these minor merger events. Nonetheless, it is encouraging 
that the shapes and the peaks of the SFR and BHAR densi- 
ties we find are significantly different, clearly indicating that 
the BHAR is not simply following the evolution of the SFR 
with time. Finally, by comparing the simulations with and 
without AGN feedback we conclude that the SFR density 
starts to be noticeably reduced by BH heating for z < 4, 
and at 2 = 1, this reduction is of order ~ 20%. 

In Di Matteo et al. (2007), similar results for the BH 
mass density and BHAR density are obtained for a cosmo- 
logical box performed at somewhat higher numerical resolu- 
tion, further supporting the trends we find here. The peak 
of the BHAR density appears moderately broader here than 
in Di Matteo et al. (2007), which could be a result of our 
smaller simulation box, or simply due to different choices of 
cosmological parameters, or due to the absence of galactic 
winds in Figure 12. As we shall discuss in detail in Section 6, 
the 'radio mode' does not play a significant role at high red- 
shift and hence is unlikely to modify the shape of the BHAR 
evolution at the peak of the quasar activity. 

In the left-hand panel of Figure 13 we show how the 
BHAR evolves with redshift as a function of the BH mass. 
We have already considered an analogous plot for our galaxy 
cluster simulations in Figure 7. While the two plots are con- 
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sistent for the overlapping range of BH masses, the most 
striking difference can be seen at the high BH mass end; in 
the galaxy cluster simulation we saw an upturn of the BHAR 
at low redshifts for very massive BHs, while here instead the 
BHAR keeps decreasing with increasing BH mass at all red- 
shifts considered. This can be understood as a consequence 
of the limited size of the cosmological box. As a result of the 
small volume, rare big objects are absent, and hence massive 
BHs residing in cluster-sized objects are missing as well. 

The characteristic time of BH growth, as a function 
of the BH mass, is illustrated in the right-hand panel of 
Figure 13. For each BH, we have defined the characteristic 
growth time as the ratio of its mass to its accretion rate 
at the given redshift, normalized to the Hubble time at the 
present epoch. It can be seen that the BH growth time is 
shortest at high redshifts and for the smallest BHs. At z = 1, 
BHs with mass of order of W^Mq reach a characteristic 
growth time of the order of the Hubble time for the first 
time. At different redshifts, the log tc — log Mbh curve has 
a similar shape, which is consistent with results obtained for 
the local galax;y population (Heckman et al., 2004). 

5.2 EfFects on the intergalactic medium 

In Figure 14, we show projected mass- weighted temperature 
and metallicity maps of our cosmological box at ^ = 1. In 
the upper row, the run with cooling and star formation only 
is plotted, while the other two rows give results for simula- 
tions where the AGN feedback was included. In the lower 
row, feedback by galactic winds was additionally included. 
The positions of BHs more massive than 2 x 10'^/i~^M© 
are marked with black dots. It can be seen that galactic 
winds modify the gas temperature and metallicity properties 
more significantly than the BH feedback alone. As expected, 
galactic winds are especially important in low mass systems, 
given that they can expel gas from objects with small escape 
speeds relatively easily due to their shallow potential wells. 
Galactic winds can therefore pollute a sizable fraction of the 
volume of the intergalactic medium (IGM) with metals, to 
a level of 10~^ — 10~^ of the Solar metallicity, in agreement 
with a number of observational findings (e.g. Ranch et al., 
1997; Cowie & Songaila, 1998; Schaye et al., 2000; Ellison 
et al., 2000). 

However, as can be seen from Figure 14, AGN con- 
tribute to the process of metal pollution of the low den- 
sity gas, albeit with lower efficiency. If only AGN feedback 
is considered, we find that metals can be transported up 
to a distance of ~ 1.5/i~^Mpc from a galactic centre, but 
the large-scale metallicity distribution remains extremely 
patchy. Thus, the main difference relative to enrichment by 
galactic winds lies not only in the level of IGM metal enrich- 
ment, but in the volume filling factor of the metal-enriched 
regions. Nevertheless, by redshift z = S, AGN feedback is 
making an important contribution to the enrichment of the 
IGM, where the characteristic size of high metallicity regions 
is of the order of ~ 1 hr^Mpc. However, we should point out 
that the AGN impact on the metal distribution found in our 
simulations is probably a lower limit, due to the fact that we 
cannot account for the contribution of very bright quasars 
powered by massive BHs at high redshifts. 

Interestingly, even though AGN heating is affecting the 
gas metallicity of the IGM, it hardly aifects its mean tem- 



perature, which is probably a refiection of the small volume 
filling factor of AGN-driven outflows. In contrast, galactic 
winds raise the intergalactic temperature significantly. Since 
this can delay the formation of small galaxies, the number 

and the mass of galaxies at a given epoch is reduced as well, 
which in turn slows the growth of the BH mass density. 

5.3 Galaxy properties and evolution 

We now study the question how the properties of galaxies are 
affected by BH feedback, and how the BH masses relate to 
their hosts. For this purpose, we identify the galaxies present 
in our cosmological box at different redshifts using a special 
group finder, and compute a number of properties for them, 
e.g. their stellar mass M«, stellar velocity dispersion cr«, their 
total SFR, and their colours in the rest-frame SDSS bands. 
This also allows us to look for correlations with their central 
BH masses. Our scheme for identifying simulated galaxies as 
isolated groups of stars is based on a simplified variant of the 
SUBFIND algorithm (Springel et al., 2001a). We first com- 
pute an adaptively smoothed baryonic density field for all 
stars and gas particles, allowing us to robustly identify cen- 
tres of individual galaxies as isolated density peaks. Using 
SUBFIND's approach for growing these centres by attaching 
nearby particles of ever lower density, we then obtain a list 
of cleanly separated galaxies. By restricting the set of gas 
particles processed in this way to those that have at least a 
density of 0.01 times the threshold density for star formar 
tion, we avoid that spatially separated galax:ies axe linked 
together, while all the stars present in the simulation are 
assigned to the different galaxies. A gravitational unbinding 
procedure is then not necessary in this procedure. Note that 
we are not interested in the gaseous components of galaxies 
here. We only include the gas particles in the galaxy finding 
procedure because the dense, cold gas that most galaxies 
contain makes the method more robust. 

In Figure 15, we show the relation between BH mass and 
stellar velocity dispersion a» (left-hand panels), and between 
BH mass and M, (right-hand panels) , at three different red- 
shifts, z = 1, 2 and 3. Here we evaluate M« and a, within 
the effective radius. Re, chosen as the half-mass radius. 
Green star symbols are for the run without galactic winds, 
while red diamonds are for the simulation that also includes 
galactic winds. The results shown are from our intermedi- 
ate resolution box, but we have verified good convergence 
in these quantities by comparing with the R3 box. With 
the dashed lines we overplot the locally observed Mbh — c* 
and Mbh — M» relations, as determined by Tremaine et al. 
(2002) and Haring & Rix (2004), respectively. A number of 
interesting features are noteworthy in this figure: 

(i) Both relations show some evolution with redshift that 
however seems most prominent at the low BH mass end and 
for the run with galactic winds. The most massive objects 
found at every epoch stay very close to the local relation 
over the redshift interval considered. 

(ii) At high redshifts, galactic winds give rise to less mas- 
sive galaxies with somewhat lower cr,, and for a given M* 
or (T*, BHs are less massive when the winds are "switched- 
on". However, by redshift z = 1 and for Mbh > 5 x 10*^ M©, 
the central BH masses are comparable for a given stellar 
mass of a galaxy, regardless of the presence of winds. This 



© 0000 RAS, MNRAS 000, 000-000 



20 Sijacki et al. 





X [ /i"' kpc ] X [ /i"' kpc ] 

Figure 14. Projected mass-weighted temperature maps (left-hand panels) and gas metallicity maps (right-hand panels) at z = 1 of the 
R2 run. Top row; simulation with cooling and star formation only; middle row: simulation with additional BH model; bottom row: run 
with BH model and galactic winds "switched-on" . The positions of BH particles more massive than 2 X WH-'^Mq are marked with 
black dots. 
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Figure 15. BH mass-stellar velocity dispersion relation (left-hand panels) and BH mass— stellar mass relation (right-hand panels), at 
redshifts z = 1, 2 and 3. Green star symbols denote the run without galactic winds, while red diamonds are for a run where galactic 
winds have been included as well. Both simulations have been performed at the resolution of the R2 cosmological box. The dashed lines 
give the locally observed relationships between the considered quantities, as determined by Tremaine et al. (2002) and Haring & Rix 
(2004). 
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Figure 16. Specific SFR (left-hand panel) and g — r colour (right-hand panel) for simulated galaxies as a function of stellar mass at 
2 = 1. Blue dots and continuous contours are for the run with cooling and star formation only. Green dots and long-dashed contours 
indicate the case where our BH model was included, while red dots and short-dashed contours are for the simulation which in addition 
contains galactic winds. 



suggests that the feedback by galactic winds is delaying BH 
growth by expelling significant amounts of gas at high red- 
shift. However, this gas becomes available for BH accretion 
at latter epochs, when it has been reincorporated from the 
IGM back into the galaxies. 

(iii) For BH masses smaller than 5 x 10® Mq (as indicated 
with the arrow in Figure 15) there is a tilt in both relations, 
especially noticeable for the Mbh — M, relation. This tilt 
is significantly reduced in the run with galactic winds. This 
can be understood by considering that the AGN feedback 
in these small galaxies is very modest, and does not modify 
the properties of the hosts significantly. On the other hand, 
galactic winds are very efficient for these low mass objects, 
reducing their stellar mass and bringing them into better 
agreement with the expected relation. However, we caution 
that this result of our modelling may be sensitive to our pre- 
scription for BH seeding and the early initial growth, which 
could be affected by numerical resolution effects. If small 
mass galaxies could grow somewhat bigger BHs in their cen- 
tre, they would also be more affected by BH feedback so that 
it may not be necessary to invoke galactic winds to match 
the low mass end of the Mbh — o" relation. In any case, it is 
interesting that galactic winds shift these smaller objects in 
the desired direction while at the same time more massive 
galaxies are unaffected and still lie on the local relation. 

(iv) Finally, given that the simulated A/bh — M* and 
Mbh — o", relations at z = 1 match the locally observed 
relations well it would be interesting to explore how the 
model will evolve further until z — 0. It seems quite hkely, 
especially for more massive galaxies, that not much evolu- 
tion will occur for z < 1, given that these galaxies have 
a reduced star formation and depleted gas content due to 
the AGN feedback. However, considering the complex inter- 
play between BHs and galactic winds, a reliable answer of 
this question will require explicit numerical simulations of a 



bigger cosmological box at comparable or better resolution, 
which we will tackle in future work. 

In Di Matteo et al. (2007), a measurement of the 
Mbh — M, and A^bh ~ o"* relations is presented as well, 
for an overlapping redshift range. Their results extend to 
somewhat higher BH mass, due to their larger box-size and 
their higher normalization as of the power spectrum. In gen- 
eral, we find a reassuring agreement with the findings of Di 
Matteo et al. (2007), indicating that the Mbh — M, and 
Mbh — o", relations are not significantly affected by the bub- 
ble feedback or the underlying cosmology, as we will discuss 
latter on in more detail. Only for the Mbh — a* relations 
at z = 3, our measurements lie somewhat below those of Di 
Matteo et al. (2007), a trend that we attribute to the fact 
that we required dark matter halos to be five times more 
massive than Di Matteo et al. (2007) before we endowed 
them with a seed black hole, which may have delayed their 
early growth. 

In Figure 16, we show the specific SFRs and the g ~ r 
colours of our simulated galaxies at 2: = 1 as a function 
of their stellar mass. Blue dots and blue continuous con- 
tours show the results when only cooling and star forma- 
tion is considered. Green dots and long-dashed contours give 
the results for the run with AGN feedback, while red dots 
and short-dashed contours represent the case with additional 
galactic winds. While AGN feedback does not affect the bulk 
of the population significantly, it is evident that the stellar 
masses and the SFRs of most massive galaxies are reduced. 
When feedback by galactic winds is included as well this 
trend is even more pronounced, and extends to the lower 
mass systems. Note also that the galaxy colours change in 
our BH model - a number of galaxies become much red- 
der, forming a red sequence, as it is shown in the right-hand 
plot of Figure 16. However, if the galactic wind model is in- 
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Figure 17. Distribution of the BHAR, expressed in MQ/yr, as a function of accretion rate in Eddington units, for three different redshifts, 
as indicated on the panels. The left-hand panel shows our measurements for the cosmological box simulation at the intermediate resolution 
{R2), while the right-hand panel gives analogous results for the g676 galaxy cluster simulation. The vertical dashed line denotes the 
transition adopted in our model between BHs accreting in the "quasar" mode and those accreting in the radiatively inefficient "radio" 
mode. 



eluded, the clearly defined red sequence disappears. Instead, 
all galaxies are in general less massive and redder. The fact 
that the red sequence becomes largely unpopulated in the 
run with galactic winds is not necessarily problematic, how- 
ever. A larger cosmological box may again fill this region of 
the diagram with a population of more massive galaxies. 

Finally, in order to assess whether the colours of our 
galaxies in the AGN feedback model are realistic, we have 
computed the mean u — r colour of the red sequence at z — 1, 
obtaining u — r ~ 2.2. We can compare this result to the 
work of Bell et al. (2004) on the COMBO-17 survey. Using 
the conversion between U — V and u ~ r colours suggested 
by these authors, we find that the colour of our simulated 
galaxies along the red sequence is in good agreement with 
the early- type galaxies from the COMBO-17 survey. 

5.4 Dependence on cosmological parameters 

In order to evaluate how sensitively our BH model de- 
pends on the underlying cosmological model, we have per- 
formed two additional simulations at intermediate resolution 
(R2) with exactly identical BH parameters. In one case we 
have adopted cosmological parameters consistent with the 
WMAP Ist-year data analysis, while for the other run we 
have selected the updated parameters of WMAP's 3rd year 
data release. The main difference between the corresponding 
cosmological models lies in a reduction of the normalization 
parameter erg (0.9 0.75), a lowering of the matter density 
Q,m (0.3 0.26) and in the introduction of a tilt in the 
primordial power spectrum slope Us (1.0 — > 0.938). 

Our results show that BHs form earlier in the WMAP-1 
cosmology, as expected from the fact that host halos above a 
given mass threshold appear earlier in this model due to the 
higher normalization, and are therefore seeded with a BH 



at higher z. Also, the number density of BHs is higher and 
there are more high mass BHs at lower redshift with respect 
to the WMAP-3 run. Consequently, the comoving BH mass 
density is somewhat higher in the WMAP-1 cosmology. At 
2 = 1, it is 3.9 X lO'^MoMpc"^, while for WMAP-3 we ob- 
tain 2.2 X IO^MqMpc""^ at the same numerical resolution. 
The BHAR density is higher in the WMAP-1 model as well, 
with its peak shifted to a slightly earlier redshift of ~ 3.5. 
However, for lower redshifts this is compensated by a steeper 
decline, such that at z = 1 the BHAR densities have a very 
similar value in both cosmologies. These trends can simply 
be understood as a result of "delayed" structure formation in 
the WMAP-3 cosmology. Increasing erg will let the BHs form 
earlier, evolve faster, and reach their peak activity at higher 
redshift. Interestingly, within the scatter, the Mbh — o", and 
Mbh — Mt relations are however not significantly affected 
by the change in the cosmological parameters. 



6 RELEVANCE OF RADIO VERSUS QUASAR 
FEEDBACK MODE 

In this section, we explore the relative importance of our two 
different modes of black hole growth and quasar feedback for 
the build-up of the cosmological black hole mass density. For 
this purpose, we show in Figure 17 the distribution function 
of the net black hole mass accretion rate in M0/yr, as a 
function of the Eddington-normalized accretion rate. In the 
left-hand panel, we give results for our cosmological box R2 
at redshifts z; = 1, 2 and 3, while on the right-hand side 
we show results for the g676 galaxy cluster simulation at 
different redshifts. 

It can be seen that at z = 3 the bulk of the BH growth 
occurs in the "quasar" mode, where the accretion rate is 
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more than 0.01 times the Eddingtoii rate. In fact, the "ra- 
dio" mode contributes less than 2% to the total BHAR at 
this epoch, and for the g676 cluster simulation, the contri- 
bution from growth in this radiatively inefficient mode is 
similarly small (~ 1%). At lower rodshifts, BHs accreting 
at low rates make up for an increasingly higher fraction of 
the total BHAR density (which however declines with time). 
For the R2 simulation at z = 2 they contribute ~ 13%, and 
at « = 1 they are at ~ 37%. Instead, for the g676 galaxy 
cluster, these numbers arc somewhat higher, with 65% at 
z = 1 and 96% at 2; = coming from BHs in a low ac- 
cretion state. Thus, while in the cosmological box BHs in 
the quasax regime are the main channel of BH growth at all 
epochs considered, for the galaxy cluster run there appears 
to be a transition epoch at 2; ^ 1, below which BHs in the 
"radio" mode are driving the BH growth at low redshifts. 
This difference between the BHAR distributions for the R2 
run and the g676 cluster simulation can be explained by the 
presence of a dominating, very massive BH in the centre 
of the main progenitor of the cluster simulation, while the 
cosmological box contains a much fairer sample of the black 
hole mass function. 

However, integrating the BHAR over cosmic time for all 
BHs in the simulated volumes, we find that the bulk of BH 
mass is always grown during phases of quasar activity, i.e. in 
the high accretion rate regime. In fact, BH accreting in the 
"radio" mode contribute less than ~ 5% to the integrated 
black hole mass density in the cosmological box, while for the 
g676 cluster simulation this immbcr is ~ 20%. Performing 
the same calculations for our gl galaxy cluster simulation we 
obtain a very similar result 15%), further strengthening 
this point. Thus, we conclude that most of the BH mass is 
assembled during periods of rapid growth in the radiatively 
efficient "quasar" mode. This implies that our models are 
consistent with the Soltan argument (Soltan, 1982; Yu & 
Tremaine, 2002; Marconi et al., 2004; Merloni, 2004; Hopkins 
et al., 2006b), and shows that merger-driven quasar activity 
is still dominating the black hole assembly in our unified 
model for AGN feedback. 

Furthermore, by repeating our cosmological box simu- 
lations without "radio" mode feedback we have explicitly 
checked that the inclusion of bubble feedback has a negli- 
gible influence on the Mbh — Mt and Mbh — cr* relations. 
Only aX z < 2, the "radio" mode makes BHs slightly less 
massive due to the somewhat more efficient heating by the 
bubbles. Thus, quasar activity triggered by the merging of 
host galaxies is the primary process responsible for establish- 
ing the scaling relationships between the central supermas- 
sive BHs and their host galaxies, and this process remains 
intact and is essentially unaffected if an additional "radio" 
mode feedback is invoked in the low accretion state of BHs. 

On the other hand, our simulations also show that the 
"radio" mode feedback in the form of AGN-driven bub- 
bles is essential in massive groups and clusters. In order 
to demonstrate this point directly, we have rerun our g676 
galaxy cluster simulation by switching off the radio channel 
of feedback mode, only allowing the ordinary quasar activ- 
ity to proceed until z = 0. This results in simulated ICM 
properties that are in disagreement with observations. In 
particular, the gas temperature profile keeps rising towards 
the very centre. Thus the intermittent nature of our bubble 
feedback appears necessary in this respect, as well of course 



for explaining the observed phenomena of radio galaxies and 
X-ray cavities in clusters of galaxies. 



7 DISCUSSION AND CONCLUSIONS 

In this study, we have proposed a new model for follow- 
ing BH growth and feedback in cosmological simulations of 
structure formation. Within our prescription, BH seeds are 
introduced at early cosmic times in the centres of all halos 
once their mass exceeds a certain threshold value for the first 
time. The BH seeding is accomplished by frequently running 
a FoF algorithm on the fly in our simulation code, which 
identifies newly formed halos and their properties. The BHs 
themselves are represented as coUisionless sink particles that 
can grow their mass via two channels, by (1) gas accretion 
with an accretion rate estimated by a simple Bondi formulae, 
and (2) via mergers with other BHs that arc close enough. 
Motivated by the observed phenomenology of radio galax- 
ies and quasars, and by analogy with X-ray binaries, we 
have assumed that the AGN feedback that accompanies gas 
accretion can be decomposed into two physically distinct 
modes: BHs accreting at high accretion rates in terms of 
their Eddington limit are operating in a radiatively efficient 
"quasar" regime, while BHs accreting at much lower rates 
live in a radiatively inefficient "radio" mode, which however 
produces significant mechanical luminosity that gives rise to 
jet-inflated bubbles. 

In our initial tests of the model, we have applied it to 
simulations of isolated galaxy clusters consisting of a static 
dark matter potential and gas in hydrostatic equilibrium. 
Considering a range of cluster masses from lO"/i"^M0 to 
10^''h^^MQ, we have found that the model leads to a self- 
regulated AGN feedback that produces realistic ICM prop- 
erties and prevents the overcooling problem in the central 
cluster regions. At the same time, it yields reasonable BH 
masses and accretion rates. 

We have then applied our model to fully self-consistent 
cosmological simulations of galaxy cluster formation, where 
the seeding of BHs and their subsequent growth is followed 
self-consistently. Together with our work in Di Matteo et al. 
(2007), these simulations are the first cosmological hydro- 
dynamical models that simultaneously follow the growth of 
structure in the dark matter as well as the baryonic physics 
of star formation and AGN feedback. These simulations pro- 
vide information on how the BHs are spatially distributed 
in the simulated volume, how they grow with cosmic time 
as a function of their environment, and how they influence 
the surrounding medium and their host galaxies. 

Our analysis of these calculations has focused on the 
history of the BH in the central cD galaxy. We have shown 
that while essentially all the mass in supermassive BHs orig- 
inates in gas accretion, a large fraction of the mass of the BH 
in the cD grows from mergers with other BHs, as opposed to 
being all accreted by a single massive progenitor. This sug- 
gests that mergers of supermassive BHs arc an important 
factor for building up very massive BHs. In line with our 
findings for the isolated cluster simulations we found that 
feedback from BHs is preventing the formation of excessive 
cooling flows out of hot cluster atmospheres, and is chang- 
ing the central ICM properties substantially in the process. 
These changes are all acting to bring the simulated temper- 
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aturc and entropy profiles into mucli better agreement witli 
observational findings of cool core clusters (e.g. De Grandi &: 
Molendi, 2002; Pratt & Arnaud, 2003; Ponman at al., 2003; 
Vikhlinin at al., 2005; Bauar at al., 2005; Dunn & Fabian, 
2006; Sanderson ct al., 2006; Pratt et al., 2007). Moreover, 
the stellar properties of the central galaxy are affected as 
wall, with the forming cD galaxies being less massive, and 
having older and redder stellar populations which are con- 
sistent with the results found in recent studies of a large 
sample of BCGs in the SDSS survey (von der Linden et al., 
2007). 

We have also studied our unified AGN feedback model 
in simulations of homogeneously sampled cosmological vol- 
umes, thereby investigating how successful it is on the scale 
of galaxies, and with respect to the mean properties of the 
cosmic population of BHs. In our simulation box of size 
25 h~^Mpc on a side, we were able to resolve galaxies with 
stellar masses greater than ~ lO^M©. We have found that 
our model produces a comoving BH mass density equal to 
2.8 X lO''M0Mpc^'', whicii is in very good agreement with 
observational estimates (Fabian & Iwasawa, 1999; Merritt & 
Ferrarese, 2001; Yu & Tremaine, 2002; Cowie et al., 2003). 
Moreover, the AGN feedback influences the properties and 
the formation of the host galaxies and vice versa, and this 
mutual coupling establishes Mbh — c* and Mbh — Af« relar 
tionships already at early epochs. Again, these relations are 
in broad agreement with observations, at least for the more 
massive and better resolved systems. Interestingly, a model 
that besides BH feedback includes feedback from supernova- 
driven galactic winds produces a better match to the ex- 
pected relationships at the low mass end. At the same time, 
this model also leads to a much more widespread enrichment 
of the IGM, which appears required by the data on quasar 
absorption line systems. Additionally, AGN heating reduces 
the stellar mass of the most luminous galaxies identified in 
the simulated volume, and it affects their SFR and colours as 
well. The galaxies in simulations with AGN feedback show 
a clearly defined red sequence already at z = 1. 

Finally, we have explored the relative importance of "ra- 
dio" versus "quasar" mode. We have found that the bulk of 
the BH growth occurs at high accretion rates, corresponding 
to radiatively efficient AGN activity. While the relative im- 
portance of the "radio" mode grows towards late times, and 
becomes large in clusters of galaxies aX z < 1, this mode of 
accretion is unimportant for the total black hole mass den- 
sity today. This also implies that our model is consistent 
with the Soltan argument. We also found that the feedback 
of the quasar regime is the key factor for establishing the 
Mbh — cr« and Mbh — M* relationships, a process that ap- 
pears unaffected by the bubble feedback. On the other hand, 
the mechanical luminosity of AGN-driven bubbles in groups 
and cluster of galaxies is necessary in order to successfully 
reproduce the observed ICM properties. 

These results represent highly encouraging successes for 
hydrodynamical models of hierarchical galaxy formation in 
ACDM cosmologies, and provide support for the theoreti- 
cal conjecture that galaxy formation and supermassive BH 
growth are intimately linked. In spite of these successes we 
need to emphasize that our BH model clearly represents a 
drastically simplified picture of BH physics. In part, the sim- 
plifications arc driven by numerical limitations, because even 
state-of-the-art simulations of galaxy formation in cosmolog- 



ically representative regions of the Universe cannot directly 
resolve the gravitational sphere of influence of individual 
BHs. We are therefore forced to adopt a subresolution ap- 
proach for the representation of BH accretion, and also for 
the modelling of bubble feedback where the initial stages 
of bubble inflation by a jet are not resolved. However, this 
subresolution approach can still capture the essential parts 
of the physics relevant for the coupling of the BHs to their 
larger-scale environment, and this is what allows us to study 
their cosmological signiflcance for galaxy formation. 

While we have tried in this study to shed some light 
onto the entangled histories of galaxies and their central 
BHs, a number of important questions remain open for fu- 
ture work, and our model appears suitable to help answering 
them. They include: (i) The early growth of BHs - can the 
observed 'downsizing' of BH growth be reproduced by our 
model? (ii) Do the galaxy cluster scaling relations change as 
a result of the AGN feedback? (iii) Can a galaxy luminos- 
ity function with a bright end consistent with observational 
constraints be reproduced in cosmological simulations, and 
what is the respective AGN luminosity function? (iv) Do the 
quasar clustering properties come out right? 

Indeed, AGN appear to be a key ingredient of struc- 
ture formation, and future numerical simulations should be 
a very helpful theoretical tool to unveil the intriguing and 
complex picture of their interactions and growth. 
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